//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include "ephi.hpp"
#include "threading.hpp"

#define CURRENT_START 60000.0
#define CURRENT_STEP 2000.0
#define ENERGY 5000.0
#define SAMPLECOUNT 10000

struct simu : public Task
{
  simu (Statics& statics, const vect3d& pos, const vect3d& dir, prec_t energy, prec_t outmag) :
      dyn(statics, true), outmag(outmag)
  {
    dyn.setDTM(0.1, 10);
    index = dyn.inject_electron (pos, dir.normal(), energy);
  }

  void execute ()
  {
    while (1)
      {
        for ( int i = 0; i < 100; i++ )
            dyn.step();
        vect3d p = dyn.get_pos(index);
        prec_t pmag = p.magnitude();
        if ( pmag > outmag )
          {
            vect3d v = dyn.get_v(index);
            if ( (p + v.normal()*0.01).magnitude() > pmag )
                break;
          }
      }
  }

  ElectroDynamics dyn;
  prec_t outmag;
  size_t index;
};

int main (int argc, char *argv[])
{
  Ephi::setDebugLevel (Ephi::DEBUG_CRITICAL);
  prec_t current_start = CURRENT_START;
  prec_t current_step = CURRENT_STEP;
  prec_t energy = ENERGY;
  int samplecount = SAMPLECOUNT;

  prec_t radius = 0.15;
  prec_t spacing = 0.07;
  vect3d pos(0.0, 0.0, 0.0);

  prec_t realR = radius + (spacing) / PREC_SQRT2;

  prec_t outmag = radius * 3;
  outmag *= outmag;

  TaskManager tmgr(4);

  prec_t v = ElectroDynamics::eV2ms (energy, ELECTRON_MASS);

  FILE *fp = fopen ("polywell_confine.csv", "w");
  if (fp)
    {
      fprintf (fp, "RADIUS[m],SPACING[m],CURRENT[A],ENERGY[eV],SPEED[km/s],AVG(TIME)[s],AVG(PATHLEN)[m],MAX(PATHLEN)[m]\n");
      fflush (fp);
    }

  for ( int currentN = 0; currentN <= 50; currentN++ )
    {
      prec_t realI = current_start + currentN * current_step;

      prec_t tsum = 0;
      prec_t tmax = 0;

      Statics statics;
      make_polywell_cube (statics, radius, 0, spacing, realI);

      for ( int i = 0; i < samplecount; i++ )
        {
          vect3d target(prec_t_drand() * 2 * realR - realR, prec_t_drand() * 2 * realR - realR, realR);
          vect3d dir = target - pos;

          tmgr.addTask (new simu(statics, pos, dir, energy, outmag));
        }

      simu *s;
      int c = samplecount;
      while ( (s = (simu *)tmgr.run()) )
        {
          prec_t t = s->dyn.get_elapsed_time();
          tsum += t;
          if ( tmax < t )
              tmax = t;
          c--;
          fprintf (stderr, "\r%d ", c);
        }
      fprintf (stderr, " \r");

      printf ("%f,%f,%f,%f\n", prec2double(realI), prec2double(tsum / samplecount), prec2double((tsum / samplecount) * v), prec2double(tmax * v));
      if ( fp )
        {
          fprintf (fp, "%f,%f,%f,%f,%.1f,%e,%f,%f\n",
              prec2double(radius), prec2double(spacing), prec2double(realI), prec2double(energy), prec2double(v / 1000.0), prec2double(tsum/samplecount), prec2double((tsum / samplecount) * v), prec2double(tmax * v));
          fflush (fp);
        }
    }

  if (fp)
      fclose (fp);
}

